No convergence issues, largest
| mean | sd | Rhat | n.eff | |
|---|---|---|---|---|
| B.hat[1] | 149.3894907 | 1.7142280 | 1.006563 | 380 |
| B.hat[2] | 6.5280123 | 0.3644497 | 1.005039 | 500 |
| sigma.b0 | 8.6347718 | 1.3826121 | 1.004578 | 550 |
| sigma.b1 | 1.8019277 | 0.2887957 | 1.007852 | 340 |
| rho | 0.6077648 | 0.1406052 | 1.000374 | 2000 |
| sigma.y | 0.6623925 | 0.0339128 | 1.001047 | 2000 |
They agree well with lmer estimates:
| Std.Dev | Correlation | |
|---|---|---|
| Intercept | 8.081 | NA |
| Slope | 1.681 | 0.641 |
| Residual | 0.660 | NA |
| Estimate | Std. Error | t value | |
|---|---|---|---|
| (Intercept) | 149.371753 | 1.5854173 | 94.21605 |
| age | 6.525469 | 0.3363003 | 19.40370 |
I removed xi.a and got better convergence with max
| mean | sd | Rhat | n.eff | |
|---|---|---|---|---|
| mu.a | 149.3060324 | 1.5523646 | 1.000601 | 5000 |
| mu.b | 6.5169483 | 0.3434120 | 1.001047 | 5000 |
| sigma.a | 7.9999186 | 1.1396635 | 1.000886 | 5000 |
| sigma.b | 1.7210504 | 0.2592583 | 1.016492 | 160 |
| rho | 0.6134800 | 0.1271427 | 1.010672 | 260 |
| sigma.y | 0.6654763 | 0.0347933 | 1.001397 | 3100 |
| _Again w | e have good ag | reement in c | oefficient | and SD estimates._ |
owls <- read.table("http://www.math.montana.edu/~jimrc/classes/stat506/data/Owls.txt", head=TRUE)
owls$cTime <- owls$ArrvTime - 24
require(lme4)
require(splines)
## Loading required package: splines
owl.fit3 <- glmer(Ncalls ~ offset(log(BroodSize)) + FoodTrt*SexParent +
SexParent*bs(cTime, df=5) +(1 + cTime|Nest), data=owls,
family=poisson)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00271455
## (tol = 0.001, component 13)
require(R2jags, quietly=TRUE)
cat("model{
for (i in 1:n){
y[i] ~ dpois (lambda[i])
log(lambda[i]) <- offset[i] + X[i,] %*% beta[] + epsilon[i]
epsilon[i] ~ dnorm (0, tau.epsilon)
}
tau.epsilon <- pow(sigma.epsilon, -2)
sigma.epsilon ~ dunif (0, 100)
for (col in 1:K){
beta[col] ~ dnorm (0, .0001)
}
}
", file = "owlsFE.jags")
logBsize <- log(owls$BroodSize)
X <- with( owls, model.matrix(~ SexParent *(FoodTrt + bs(cTime, df = 5) )))
## dim(X) 599 * 14
owlsData <- with(owls, list("y" = Ncalls, "X" = X, "offset" = logBsize, "n" = nrow(X),"K" = ncol(X) ) ) #, "J" = length(levels(Nest))
owl.params = c("beta", "sigma.epsilon")
owl.jags1 <- jags( data = owlsData, inits = NULL, parameters.to.save = owl.params,
model.file = "owlsFE.jags", n.chains=4, n.iter=1000, progress.bar="none")
## module glm loaded
Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 12504
Initializing model
# owl.jags1
owl.fit1 <- glm(Ncalls ~ offset(log(BroodSize)) + SexParent*(FoodTrt + bs(cTime,5)), data = owls, family=quasipoisson)
glmCoef <- coef(owl.fit1)
owl.jags1Coef<- apply(owl.jags1$BUGSoutput$sims.matrix, 2, mean)[1:ncol(X)]
owl.jags1SECoef<- apply(owl.jags1$BUGSoutput$sims.matrix, 2, sd)[1:ncol(X)]
se.lmCoef <- summary(owl.fit1)$coef[,2]
coefTable <- rbind( owl.jags1Coef, glmCoef, owl.jags1SECoef, se.lmCoef)
colnames(coefTable) <- paste("b",1:14,sep="")
knitr::kable(coefTable,digits=3)
| b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 | b10 | b11 | b12 | b13 | b14 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| owl.jags1Coef | 0.724 | -1.928 | -0.950 | 0.885 | -1.660 | -0.081 | -1.113 | -1.164 | 0.069 | 2.618 | 1.796 | 3.530 | -0.444 | 3.837 |
| glmCoef | 0.618 | -1.045 | -0.545 | 1.088 | -0.785 | 0.192 | -0.380 | -0.565 | 0.050 | 1.408 | 0.827 | 2.255 | -0.830 | 2.162 |
| owl.jags1SECoef | 0.574 | 1.042 | 0.174 | 0.890 | 0.594 | 0.916 | 0.817 | 0.939 | 0.230 | 1.489 | 0.995 | 1.456 | 1.271 | 1.432 |
| se.lmCoef | 0.405 | 0.729 | 0.130 | 0.642 | 0.422 | 0.706 | 0.628 | 0.727 | 0.166 | 1.051 | 0.688 | 1.071 | 0.937 | 1.011 |
owl1.mcmc <- as.mcmc(owl.jags1)
## colnames(owl1.mcmc[[1]])
#for(ndx in 1:length(owl1.mcmc)){
## owl1.mcmc[[ndx]] <- owl1.mcmc[[ndx]][,c(1, 54, 27, 54:58)]
#}
par(mfrow=c(4,4), mar = c(3,2,2,1))
traceplot(owl1.mcmc)
Having gotten that to run, and seeing fair agreement with a fixed effects glm model, I then add the random nest effects: and intercept, and a slope over centered time.
require(R2jags, quietly=TRUE)
cat("model{
for (i in 1:n){
y[i] ~ dpois (lambda[i])
log(lambda[i]) <- offset[i] + X[i,] %*% beta[] +
a.raw[nest[i]] + b.raw[nest[i]]*time[i]
}
for (col in 1:K){
beta[col] ~ dnorm (0, .0001)
}
for(j in 1:J){
a.raw[j] ~ dnorm(0, tau.a)
b.raw[j] ~ dnorm(0, tau.b)
a[j] <- a.raw[j] - mean(a.raw)
b[j] <- b.raw[j] - mean(b.raw)
}
tau.a <- pow(sigma.a, -2)
sigma.a ~ dunif (0, 100)
tau.b <- pow(sigma.b, -2)
sigma.b ~ dunif (0, 100)
}
", file = "owlsME.jags")
owlsData2 <- with(owls, list("y" = Ncalls, "X" = X, "offset" = logBsize, "n" = nrow(X),
"K" = ncol(X), "J" = length(levels(Nest)),
"nest" = as.numeric(Nest), "time" = cTime ))
owl.param2 = c("beta", "a", "b", "sigma.a","sigma.b")#,"sigma.epsilon")
owl.inits <- function (){
list ( sigma.a = runif(1)*10, sigma.b = runif(1)*10, beta = rnorm(14,0,10))
}
owl.jags2 <- jags( data = owlsData2, inits = owl.inits, parameters.to.save = owl.param2,
model.file = "owlsME.jags", n.chains=4, n.iter=1000,progress.bar ="none")
## module glm loaded
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 13806
##
## Initializing model
owl.jags2 <- autojags(owl.jags2, n.update = 5, n.iter = 5000,progress.bar ="none")
The model did converge as maximum
The plots above show jags predicted random effects on the horizontal axis against lmer predictions on the y axis. The identity line is shown in grey. One or possibly two nests are a bit above the identity line, but I think this is good agreement.
R/jags Code
Oxboys with correlation.
require(lme4, quietly=TRUE)
require(R2jags, quietly=TRUE)
data(Oxboys, package="nlme")
cat("model {
for(i in 1:n){
y[i] ~ dnorm(y.hat[i], tau.y)
y.hat[i] <- B[Subject[i], 1] + B[Subject[i], 2] * x[i]
}
tau.y <- pow(sigma.y, -2)
sigma.y ~ dunif (0, 100)
for (j in 1:J){
B[j,1:2] ~ dmnorm(B.hat[1:2], Tau.B[,])
}
B.hat[1] ~ dnorm(0, .0001)
B.hat[2] ~ dnorm(0, .0001)
sigma.b0 ~ dunif (0, 100)
sigma.b1 ~ dunif (0, 100)
rho ~ dunif(-1,1)
Tau.B <- inverse(Sigma.B)
Sigma.B[1,1] <- sigma.b0^2
Sigma.B[2,2] <- sigma.b1^2
Sigma.B[1,2] <- sigma.b0 * sigma.b1 * rho
Sigma.B[2,1] <- sigma.b0 * sigma.b1 * rho
}",
file="OB-Corr.jags")
jags2.data = with(Oxboys, list("n" = 234,"J" = 26, "y" = height, "x" = age, "Subject" = Subject))
oxboys.inits <- function (){
list ( sigma.y = runif(1)*10, sigma.b0 = runif(1)*10, sigma.b1 = runif(1)*4)
}
OB4.params <- c ("B", "B.hat", "sigma.y", "sigma.b0", "sigma.b1", "rho")
OB4.jags <- jags(data = jags2.data, inits = oxboys.inits, OB4.params,
model.file = "OB-Corr.jags", n.chains=4, n.iter=1000, progress.bar="none")
## OB4.jags
OB4.out <- print(OB4.jags)$summary[c(53:54,57:58, 56, 59), ]
## colnames(OB4.mcmc[[1]])
OB4.mcmc <- as.mcmc(OB4.jags)
keep <- match(rownames(OB4.out), colnames(OB4.mcmc[[1]]))
for(ndx in 1:length(OB4.mcmc)){
OB4.mcmc[[ndx]] <- OB4.mcmc[[ndx]][, keep]
}
par(mfrow=c(2,3), mar = c(3,2,2,1))
traceplot(OB4.mcmc)
maxRhat4 <- max(gelman.diag(as.mcmc(OB4.jags))[[1]][,1])
require(R2jags, quietly=TRUE)
knitr::kable(OB4.out[, c(1:2,8:9)])
require(lme4, quietly=TRUE)
lmer4 <- lmer(height ~ age + (1 + age|Subject), Oxboys)
## lmer Random Effects:
vc <- VarCorr(lmer4)[[1]]
sigma <- summary(lmer4)$sigma
vc <- matrix(c(attr(vc, "stddev"), sigma,NA,attr(vc, "correlation")[2,1], NA), ncol=2, nrow=3)
dimnames(vc) = list(c("Intercept","Slope","Residual"), c("Std.Dev", "Correlation"))
knitr::kable(vc, digits = 3)
## lmer Fixed Effects:
knitr::kable(summary(lmer4)$coefficients)
Oxboys with Wishart
require(R2jags, quietly=TRUE)
cat("model {
for(i in 1:n){
y[i] ~ dnorm(y.hat[i], tau.y)
y.hat[i] <- a[Subject[i]] + b[Subject[i]] * x[i]
}
tau.y <- pow(sigma.y, -2)
sigma.y ~ dunif (0, 100)
for (j in 1:J){
a[j] <- B.raw[j,1]
b[j] <- xi.b* B.raw[j,2]
B.raw[j,1:2] ~ dmnorm(B.raw.hat[j, ], Tau.B.raw[,])
B.raw.hat[j, 1] <- mu.a.raw
B.raw.hat[j, 2] <- mu.b.raw
}
mu.a <- mu.a.raw
mu.b <- xi.b * mu.b.raw
mu.a.raw ~ dnorm(0, .0001)
mu.b.raw ~ dnorm(0, .0001)
xi.b ~ dunif(0, 100)
xi.a ~ dunif(0, 100)
Tau.B.raw ~ dwish(W[ , ], 3)
Sigma.B.raw <- inverse(Tau.B.raw[,])
sigma.a <- sqrt(Sigma.B.raw[1,1])
sigma.b <- xi.b * sqrt(Sigma.B.raw[2,2])
rho <- Sigma.B.raw[1,2] / sqrt(Sigma.B.raw[1,1] * Sigma.B.raw[2,2])
}",
file="OB-Wishart.jags")
jags5.data = with(Oxboys, list("n" = 234,"J" = 26, "y" = height, "x" = age, "Subject" = Subject, "W" = diag(2)))
J <- 26
oxboys.inits <- function (){
list (B.raw = array(rnorm(2*J),c(J,2)),mu.a.raw = rnorm(1), mu.b.raw=rnorm(1), sigma.y = runif(1)*10, Tau.B.raw = rWishart(1,3,diag(2)))
}
OB5.params <- c ("a","b", "mu.a", "mu.b","rho","sigma.y", "sigma.a", "sigma.b")
OB5.jags <- jags(data = jags5.data, inits = NULL, parameters.to.save =OB5.params,
model.file = "OB-Wishart.jags", n.chains=4, n.iter=1000, progress.bar="none")
OB5.jags <- autojags(OB5.jags,n.iter=5000,Rhat = 1.05, n.update =5, progress.bar="none" )
## OB5.jags
OB5.out <- print(OB5.jags)$summary[c(54,55,57:58, 56, 59), ]
## colnames(OB4.mcmc[[1]])
OB5.mcmc <- as.mcmc(OB5.jags)
keep <- match(rownames(OB5.out), colnames(OB5.mcmc[[1]]))
for(ndx in 1:length(OB5.mcmc)){
OB5.mcmc[[ndx]] <- OB5.mcmc[[ndx]][, keep]
}
par(mfrow=c(2,3), mar = c(3,2,2,1))
traceplot(OB5.mcmc)
mean4 <- round(mean(OB4.jags$BUGSoutput$sims.list[["deviance"]]), 3)
sd4 <- round(sd(OB4.jags$BUGSoutput$sims.list[["deviance"]]), 3)
mean5 <- round(mean(OB5.jags$BUGSoutput$sims.list[["deviance"]]), 3)
sd5 <- round(sd(OB5.jags$BUGSoutput$sims.list[["deviance"]]), 3)
maxRhat5 <- round(max(gelman.diag(as.mcmc(OB5.jags))[[1]][,1]),3)
require(R2jags, quietly=TRUE)
knitr::kable(OB5.out[, c(1:2,8:9)])
owls <- read.table("http://www.math.montana.edu/~jimrc/classes/stat506/data/Owls.txt", head=TRUE)
owls$cTime <- owls$ArrvTime - 24
require(lme4)
require(splines)
owl.fit3 <- glmer(Ncalls ~ offset(log(BroodSize)) + FoodTrt*SexParent +
SexParent*bs(cTime, df=5) +(1 + cTime|Nest), data=owls,
family=poisson)
Owls fixed
require(R2jags, quietly=TRUE)
cat("model{
for (i in 1:n){
y[i] ~ dpois (lambda[i])
log(lambda[i]) <- offset[i] + X[i,] %*% beta[] + epsilon[i]
epsilon[i] ~ dnorm (0, tau.epsilon)
}
tau.epsilon <- pow(sigma.epsilon, -2)
sigma.epsilon ~ dunif (0, 100)
for (col in 1:K){
beta[col] ~ dnorm (0, .0001)
}
}
", file = "owlsFE.jags")
logBsize <- log(owls$BroodSize)
X <- with( owls, model.matrix(~ SexParent *(FoodTrt + bs(cTime, df = 5) )))
## dim(X) 599 * 14
owlsData <- with(owls, list("y" = Ncalls, "X" = X, "offset" = logBsize, "n" = nrow(X),"K" = ncol(X) ) ) #, "J" = length(levels(Nest))
owl.params = c("beta", "sigma.epsilon")
owl.jags1 <- jags( data = owlsData, inits = NULL, parameters.to.save = owl.params,
model.file = "owlsFE.jags", n.chains=4, n.iter=1000, progress.bar="none")
# owl.jags1
owl.fit1 <- glm(Ncalls ~ offset(log(BroodSize)) + SexParent*(FoodTrt + bs(cTime,5)), data = owls, family=quasipoisson)
glmCoef <- coef(owl.fit1)
owl.jags1Coef<- apply(owl.jags1$BUGSoutput$sims.matrix, 2, mean)[1:ncol(X)]
owl.jags1SECoef<- apply(owl.jags1$BUGSoutput$sims.matrix, 2, sd)[1:ncol(X)]
se.lmCoef <- summary(owl.fit1)$coef[,2]
coefTable <- rbind( owl.jags1Coef, glmCoef, owl.jags1SECoef, se.lmCoef)
colnames(coefTable) <- paste("b",1:14,sep="")
knitr::kable(coefTable,digits=3)
owl1.mcmc <- as.mcmc(owl.jags1)
## colnames(owl1.mcmc[[1]])
#for(ndx in 1:length(owl1.mcmc)){
## owl1.mcmc[[ndx]] <- owl1.mcmc[[ndx]][,c(1, 54, 27, 54:58)]
#}
par(mfrow=c(4,4), mar = c(3,2,2,1))
traceplot(owl1.mcmc)
Owls Mixed
require(R2jags, quietly=TRUE)
cat("model{
for (i in 1:n){
y[i] ~ dpois (lambda[i])
log(lambda[i]) <- offset[i] + X[i,] %*% beta[] +
a.raw[nest[i]] + b.raw[nest[i]]*time[i]
}
for (col in 1:K){
beta[col] ~ dnorm (0, .0001)
}
for(j in 1:J){
a.raw[j] ~ dnorm(0, tau.a)
b.raw[j] ~ dnorm(0, tau.b)
a[j] <- a.raw[j] - mean(a.raw)
b[j] <- b.raw[j] - mean(b.raw)
}
tau.a <- pow(sigma.a, -2)
sigma.a ~ dunif (0, 100)
tau.b <- pow(sigma.b, -2)
sigma.b ~ dunif (0, 100)
}
", file = "owlsME.jags")
owlsData2 <- with(owls, list("y" = Ncalls, "X" = X, "offset" = logBsize, "n" = nrow(X),
"K" = ncol(X), "J" = length(levels(Nest)),
"nest" = as.numeric(Nest), "time" = cTime ))
owl.param2 = c("beta", "a", "b", "sigma.a","sigma.b")#,"sigma.epsilon")
owl.inits <- function (){
list ( sigma.a = runif(1)*10, sigma.b = runif(1)*10, beta = rnorm(14,0,10))
}
owl.jags2 <- jags( data = owlsData2, inits = owl.inits, parameters.to.save = owl.param2,
model.file = "owlsME.jags", n.chains=4, n.iter=1000,progress.bar ="none")
owl.jags2 <- autojags(owl.jags2, n.update = 5, n.iter = 5000,progress.bar ="none")
owl2.mcmc <- as.mcmc(owl.jags2)
## colnames(owl2.mcmc[[1]])
for(ndx in 1:4){
owl2.mcmc[[ndx]] <- owl2.mcmc[[ndx]][,c(1, 27:28, 54:55, 68:71 )]
}
par(mfrow=c(3,3), mar = c(3,2,2,1))
traceplot(owl2.mcmc)
maxOwlRhat <- round(max(gelman.diag(as.mcmc(owl.jags2))[[1]][,1]),3)
Owl Plots
require(lme4, quietly=TRUE)
require(R2jags, quietly=TRUE)
par(mfrow=c(1,3), mar = c(2,2,2,1))
owl.jags2Coef<- apply(owl.jags2$BUGSoutput$sims.matrix[,55:68], 2, mean)[1:ncol(X)]
owl.jags2SECoef<- apply(owl.jags2$BUGSoutput$sims.matrix[,55:68], 2, sd)[1:ncol(X)]
glmerCoef <- summary(owl.fit3)$coef[,1]
se.glmerCoef <- summary(owl.fit3)$coef[,2]
plot(owl.jags2Coef, glmerCoef, main = "Fixed Effects")
abline(0,1, col = "grey")
glmerREs <- ranef(owl.fit3)[[1]]
jagREs <- summary(as.mcmc(owl.jags2))$statistics[1:54,1]
plot(jagREs[match( paste("a[",1:27,"]",sep=""),names(jagREs)[1:27])], glmerREs[,1], main = "Random Intercepts")
abline(0,1, col = "grey")
plot(jagREs[match( paste("b[",1:27,"]",sep=""),names(jagREs)[1:27+27])+27], glmerREs[,2], main = "Random Slopes")
abline(0,1, col = "grey")